function [Y, var_name, pos, L_base]=simulate_model_general(y_T,r,g_N,d,w_lag,p_lag,N_periods,dp_policy,h_policy,grid,par,Cum_prob,moment_dummy)
% Copyright (C) 2019-2023 Benjamin Born, Francesco D'Ascanio, Gernot J. Mueller, Johannes Pfeifer
%
% This is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% It is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
% 
% For a copy of the GNU General Public License,
% see <http://www.gnu.org/licenses/>.

if nargin < 13
    moment_dummy=0;
end
    
pos.r_ann=1;
var_name{pos.r_ann,1}='interest rate, $r_{ann}$';
var_name{pos.r_ann,2}='$r_{ann}$';
pos.y_T=length(fieldnames(pos))+1;
var_name{pos.y_T,1}='tradable output, $y^T$';
var_name{pos.y_T,2}='$y^T$';
pos.c=length(fieldnames(pos))+1;
var_name{pos.c,1}='final consumption, $c$';
var_name{pos.c,2}='$c$';
pos.c_T=length(fieldnames(pos))+1;
var_name{pos.c_T,1}='tradable consumption, $c^T$';
var_name{pos.c_T,2}='$c^T$';
pos.c_N=length(fieldnames(pos))+1;
var_name{pos.c_N,1}='nontradable consumption, $c^N$';
var_name{pos.c_N,2}='$c^N$';
pos.g_N=length(fieldnames(pos))+1;
var_name{pos.g_N,1}='government consumption, $g^N$';
var_name{pos.g_N,2}='$g^N$';
pos.y_N=length(fieldnames(pos))+1;
var_name{pos.y_N,1}='nontradable output, $y^N$';
var_name{pos.y_N,2}='$y^N$';
pos.unemployment=length(fieldnames(pos))+1;
var_name{pos.unemployment,1}='unemployment rate, $1-h$';
var_name{pos.unemployment,2}='$1-h$';
pos.D=length(fieldnames(pos))+1;
var_name{pos.D,1}='debt, $d$';
var_name{pos.D,2}='$d$';
pos.D_y_annual=length(fieldnames(pos))+1;
var_name{pos.D_y_annual,1}='Ann. Debt to GDP, $d/(4(y^T+py^N))$';
var_name{pos.D_y_annual,2}='$d/(4(y^T+py^N))$';
pos.w=length(fieldnames(pos))+1;
var_name{pos.w,1}='real wage, $w$';
var_name{pos.w,2}='$w$';
pos.devaluation_rate=length(fieldnames(pos))+1;
var_name{pos.devaluation_rate,1}='devaluation rate, $\epsilon$';
var_name{pos.devaluation_rate,2}='$\epsilon$';
pos.rer=length(fieldnames(pos))+1;
var_name{pos.rer,1}='real exchange rate, $RER$';
var_name{pos.rer,2}='$RER$';
pos.t_y=length(fieldnames(pos))+1;
var_name{pos.t_y,1}='tradable output to GDP ratio, $y^T/(y^T+py^N)$';
var_name{pos.t_y,2}='$y^T/(y^T+py^N)$';
pos.p=length(fieldnames(pos))+1;
var_name{pos.p,1}='relative price, $p^N$';
var_name{pos.p,2}='$p^N$';
pos.gross_wage_infl=length(fieldnames(pos))+1;
var_name{pos.gross_wage_infl,1}='nominal wage inflation, $\Delta w$';
var_name{pos.gross_wage_infl,2}='$\Delta w$';
pos.wage_infl=length(fieldnames(pos))+1;
var_name{pos.wage_infl,1}='annualized nominal wage inflation, $\Delta w$';
var_name{pos.wage_infl,2}='$(\Delta w)^4-1$';
pos.trade_balance=length(fieldnames(pos))+1;
var_name{pos.trade_balance,1}='trade balance, $y_T-c_T$';
var_name{pos.trade_balance,2}='$y_T-c_T$';
pos.y_aggregate=length(fieldnames(pos))+1;
var_name{pos.y_aggregate,1}='aggregate output, $y_T+py_N$';
var_name{pos.y_aggregate,2}='$y$';
pos.g_share=length(fieldnames(pos))+1;
var_name{pos.g_share,1}='government spending share, $g/y$';
var_name{pos.g_share,2}='$g/y$';
pos.price_infl=length(fieldnames(pos))+1;
var_name{pos.price_infl,1}='annualized nominal price inflation, $\Delta p$';
var_name{pos.price_infl,2}='$(\Delta p)^4-1$';
% pos.wcpi=length(fieldnames(pos))+1;
% var_name{pos.wcpi,1}='wage in terms of final output';
% var_name{pos.wcpi,2}='?';

Y=NaN(N_periods,length(fieldnames(pos)));

%p_lag=1;

for period_iter=1:N_periods
    % get exogenous state indices
    exo_index       = find(grid.y_level==y_T & grid.r_level==r & grid.g_level==g_N);
    %get resulting debt choice
    dp = dp_policy(exo_index,grid.d_grid==d,grid.w_grid_vec==w_lag);
    
    c_T = y_T + dp./(1+r) - d;
    %float version has h=loaded_PFI.hbar
    w_full = (1-par.omega)/par.omega*((par.hbar^par.alfa - g_N)./c_T).^(-1/par.xi) *par.alfa*par.hbar^(par.alfa-1);%full-employment real wage rate in terms of tradables
    h = h_policy(exo_index,grid.d_grid==d,grid.w_grid_vec==w_lag); %hours
    if h>=0.999
        h=1;
    end
    w = (1-par.omega)/par.omega.*((h.^par.alfa - g_N)./c_T).^(-1/par.xi)*par.alfa.*h.^(par.alfa-1);  %real wage
    % force onto grid
    w(w<grid.w_grid_vec(1)) = grid.w_grid_vec(1);
    w(w>grid.w_grid_vec(end)) = grid.w_grid_vec(end);
    w_index = round((log(w)-log(grid.w_grid_vec(1)))/grid.w_grid_step_size)+1;
    w = grid.w_grid_vec(w_index);

    y_N = h.^par.alfa; %nontradable output
    c_N = y_N - g_N; %nontradable consumption
    c = (par.omega*c_T.^(1-1/par.xi) + (1-par.omega) * c_N.^(1-1/par.xi)).^(1/(1-1/par.xi)); %total consumption
    p = (1-par.omega)/par.omega*(c_N./c_T).^(-1/par.xi); %relative price of nontradables in terms of tradables
    
    epsilon_rate = max(1,(par.gama*w_lag./w_full).^par.phi_e); %depreciation rate
    
    %inflation rates
    real_wage_inflation=w./w_lag;
    price_inflation=epsilon_rate.*((par.omega^par.xi+(1-par.omega)^par.xi*p.^(1-par.xi))./...
        (par.omega^par.xi+(1-par.omega)^par.xi*p_lag.^(1-par.xi))).^(1/(1-par.xi));
    nominal_wage_inflation=real_wage_inflation.*price_inflation;
    
    pc = 1/par.omega * (c_T./c).^(1/par.xi); 
    
    rer = (par.omega^par.xi + (1-par.omega)^par.xi*p.^(1-par.xi)).^(-1/(1-par.xi)); %real exchange rate
    t_y = (y_T)./(y_T+p.*y_N); %tradable output to output ratio
       
    Y(period_iter,pos.r_ann)=r*4;
    Y(period_iter,pos.y_T)=y_T;
    Y(period_iter,pos.c)=c;
    Y(period_iter,pos.c_T)=c_T;
    Y(period_iter,pos.c_N)=c_N;
    Y(period_iter,pos.g_N)=g_N;
    Y(period_iter,pos.y_N)=y_N;
    Y(period_iter,pos.unemployment)=1-h;
    Y(period_iter,pos.D)=dp;
    Y(period_iter,pos.y_aggregate)=y_T+p.*y_N;    
    Y(period_iter,pos.D_y_annual)=(dp./(4*Y(period_iter,pos.y_aggregate)));
    Y(period_iter,pos.g_share)=(p.*g_N)/Y(period_iter,pos.y_aggregate);    
    Y(period_iter,pos.w)=w;
    Y(period_iter,pos.devaluation_rate)=(100*(epsilon_rate.^4-1))';%JP: why annualized?
    Y(period_iter,pos.rer)=rer;
    Y(period_iter,pos.t_y)=t_y;
    Y(period_iter,pos.p)=p;
    Y(period_iter,pos.gross_wage_infl)=nominal_wage_inflation;
    Y(period_iter,pos.wage_infl)=100*(nominal_wage_inflation.^4-1);
    Y(period_iter,pos.price_infl)=100*(price_inflation.^4-1);

    Y(period_iter,pos.trade_balance)=y_T-c_T;
    
    %% get new random draw for ex. state
    rand_number=rand;
    L_base = sum(Cum_prob(exo_index,:)<rand_number)+1;
    
    y_T = grid.y_level(L_base);
    r = grid.r_level(L_base);
    g_N = grid.g_level(L_base);
     
    % update endogenous states
    p_lag = p;
    w_lag = w;
    d = dp;
    
end

if moment_dummy
    %% check targets
    fprintf(' mean debt-to-output ratio is: ');
    fprintf('%5.4f \n ',  mean(Y(:,pos.D_y_annual)));
    
    fprintf('mean tradable output share is: ');
    fprintf('%5.4f \n ',  mean(Y(:,pos.t_y)));
end